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Abstract 

Range migrations in response to climate change, invasive species and the emergence 
of novel ecosystems highlight the importance of temporal turnover in community 
composition as a fundamental part of global change in the Anthropocene. Temporal 
turnover is usually quantified using a variety of metrics initially developed to capture 
spatial change. However, temporal turnover is the consequence of unidirectional 
community dynamics resulting from processes such as population growth, colonisation 
and local extinction. Here, we develop a framework based on community dynamics, and 
propose a new temporal turnover measure. A simulation study and an analysis of an 
estuarine fish community both clearly demonstrate that our proposed turnover measure 
offers additional insights relative to spatial-context-based metrics. Our approach reveals 
whether community turnover is due to shifts in community composition or in community 
abundance, and identifies the species and/or environmental factors that are responsible 
for any change. 
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1 Introduction 


A large numbe r of m et rics have b e en de v eloped to me as ure c o mmunity dissi m ilarity 

2011 . ISsether et al. . 2013 . 


Leeendre and Gauthier. 

20 hi) in the century since 

Jaccard ( 

1901, 


method of quantifying differences in assemblage composition. This variation in the identities 
of species found in different sites, as measured by Jaccard dissimilarity and other turnover 
metrics, was dubbed /3-diversity by Whittaker (I960, 19721 ). to distinguish it from within 
assemblage diversity (known as cr-diversity). Although the emphasis, to date, has been on 
measuring /3-diversity across spatial contexts, growing concern about threat s to bi odivers ity 


u nde rlines an urgent need to quantify and understand temporal turnover ([Dornelas et al 


20131 ). Recent analyses have shown that although the a-diversity of local communities is 


not changing consistently through tim e, with many assemblages showing increas ing trends 
and others showing decreasing trends ( Vellend et ah . 20131 . Dornelas et al. . 2014), the rate 
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of change in community composition ( temporal /3-diversity) is greater than predicted by null 


models of baseline temporal turnover (IDornelas et all 12014 ISupp and Ernestl . 12014 1. 


One appr oach to quanti fyin g te mporal turnover is to use metrics such as the Jaccard and 
Bray-Curtis ( Bray and CurtisL 1957]) indices, initially developed to capture spatial change. 
However, temporal turnover has featur es, such as u nidi rectio nal change, not usually present 
in investigations of spatial /3-diversity ( Dornelas et ah . 2013h . Moreover, temporal turnover 


is the consequence of community dynamics resulting from processes such as local immigration 
and extinction, population growth and density dependence. Ideally, then, temporal turnover 
should recognise the dimensionality of temporal change and the ecological processes that lead 
to shifts in community composition through time. 

Here we formalise the concept of temporal turnover. We present a novel framework 
for measuring temporal turnover based on the dynamics of species abundance within a 
community. The turnover measure we propose is directly derived from a general form of 
a population dynamics model which describes the extent to which each population, and 
consequently the community to which these species belong, changes over time. Our framework 
highlights the fact that the concept of turnover can be dissected into two key ecological 
aspects: change in community composition and change in community size (or capacity) sensu 
its abundance (jBrown, 19811 ). In addition to providing new insights into the ecological basis 
of measures of temporal turnover, our approach has important practical applications because 
it allows the user to identify which species and/or which environmental factors play a key 
role in the change of the species community. We demonstrate how the turnover measure 
performs using a simulation and by analysing an estuarine fish community time-series from 
the Bristol Channel, UK. We also discuss the relationship between our approach and existing 
methods of a- and /3-diversities. For /3-diversity, we identify links between our measure and 
those pre viousl y deve loped to capture the nestedness and turnover components of /3-diversity 
( Baselga . 12010 . 20111201 al l. 


2 A measure of temporal turnover 

The temporal changes that all ecological communities experience are a cumulative 
consequence of changes in the abundance of each species in the community. Population 
ecology describes the dynamics of species abundance by differential equations, modelling the 
change rate of the population size of a given species. We can consider this rate of change in 
abundance as, literally, the temporal turnover of the species in question. This assumption 
is the basis of our approach. The temporal turnover of the whole community can then be 
defined as an additive effect of turnover in each species, since the total abundance of the 
community is the sum of the species abundances in the community. Keeping these points 
in mind, we propose a new measure of temporal turnover of communities, and identify two 
key ecological components of this turnover: change in community composition and change in 
total abundance. 

Throughout the paper we consider an ecological community consisting of s species whose 
abundance varies over time governed by the population’s background state at each time t. 
We use the term expected abundance to denote the background state, which is often described 
using other terms such as mean abundance, model expectation and (ideal) population size 
(a more precise definition will follow in Section [3]). We recognise that observed abundances 
reflect both expected abundance and natural variability. This recognition allows us, when we 
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have observations, to separate the effect of natural variability from the expected abundance 
in which features of turnover can be found. As such, here we construct a framework on the 
expected abundance, instead of the observations themselves. 

Let A i(t) > 0 be the expected abundance of the z-th species at time t. Note the fact that 
the observed abundance can be zero even though the expected abundance is assumed to be 
always positive. A species with very low expected abundance (close to zero) tends to be absent 
in observations most of the time. We can therefore deal with the case of species absence, 
adopting expected abundance as the basis of our framework. The expected total-abundance 
of the community is then the sum of each species’ expected abundance as A (t) = Ylt= i Xi(t). 
This fact is always valid regardless of whether the species abundances within a community 
are independent or not, and serves as the basis of our discussion below. 

We begin by examining the dynamics of the populations that make up the community. 
Consider an equation for the expected abundance of the z-th species as 


Xi(t + dt) = (1 + ai(t)dt)\i(t ), 


(1) 


where oti(t) is the instantaneous change rate at time t that drives the increase or the decrease 
of the species abundance; the abundance increases when ai(t) > 0, and decreases when 
ai(t) < 0. If the equation (Jl]) is solved, the Malthusian growth model ( Malthusl . 17981 )) 
specifies the expected abundance, Aj(f), and the instantaneous change rate, ai(t), becomes 
the growth rate of the population. The rate can also be modelled by other processes such 
as those, birth, death, immigration and emigration rates, bi(t) > 0, di(t) > 0, rrii{t) > 0 
and ei(t) > 0 respectively, like cti(t) = bi(t ) — di(t) + rrii(t) — ei(t). This of course leads 
A i(t) to another type of model. As such, in order to keep our framework below general, we 
do not specify any functional form of the expected abundance, Aj(f), at this stage, but it 
will be discussed in Section [3l Note that when those rates are offset relative to each other, 
ai(t) = bi(t ) — di(t) + rrii(t) — ei(t ) = 0, the population’s abundance is stable for the time 
increment dt, since its expected abundance becomes constant, A i(t + dt) = A,(t). 

The reasoning given above implies that change rate is a legitimate indicator to quantify, 
literally, the turnover of the z-th species, based on its abundance, for the time increment dt. 
However, as the change rate is in most cases unknown, it needs to be determined by the 
trajectory of the expected abundance, A*(f) > 0, itself as 


ai(t)dt = 


Xi(t + dt) — Xi(t) dXi(t) 


X i(t) 
= dlog (Xi(t)). 


A i{t) 


Here and elsewhere in the paper log refers to the natural logarithms. 

We now consider the community scale, represented by the expected total-abundance of 
the species in the community A (t) = applying the same analogy used for the 

population scale. Although the continuous form of the expected abundances, Aj(f),z = 
1,2is unknown, we can estimate it using observations at discrete time points. We 
therefore rewrite Equation (jT]) in a discrete form representing a community dynamics model 
over a relatively short time interval, between times t and t + h (h > 0), within which the 
difference equation is still valid. By integrating the instantaneous change rate, az(f), over the 


3 










time interval, the community dynamics can be described as 


A (t + h) = 


2=1 


rt+h 


( r 

E ( 1 + / oti(y)dv\ A i(t) 
i= 1 \ ' 

s / pt-\-h \ 

E( 1+ / dlog(Xi(v))\ Xi(t) 

i= 1 X Jt J 

Xi(t + h ) 


E ( 1 + lo S 


Aj(t) 


A i(t). 


The change rate and the instantaneous change rate may therefore hold the following 
relationship for the short time period h, 


X(t + h) — X(t) . 

A(t) 


5 Z lo s 


2 — 1 


Aj(t + h ) 

A i(t) 




( 2 ) 


where Pi(t) is the relative abundance of the i-th species at time t, that is determined as 


Pi(t) 


Xi(t) 
X (t) ' 


The collection of the relative abundances of the community — the relative abundance 
distribution p(t) = {pi(t),P2(t), ■ ■ ■ ,p s (t)} — satisfies the conditions that Pi(t) > 0 and 
£-=i Pi® = i. 

Now we can formally define the turnover measure, which we call D. From Equation (|2|), 
we define the turnover measure between times t and u, (u > t ) as 


S 

D(t : u ) = : u) 

2=1 


E lo s 


i=i 


E lo e 


2=1 


K(u) 

Xi(t) 

Pi(t ) 

Pi(u) 


Pi(t ) 


Pi(t) + log 


A(u)\ 

m) 


Di(p(t) : p(u)) + D 2 (X(t) : X(u)). 


( 3 ) 

( 4 ) 


As shown above, temporal turnover, D, is an abundance based measure, and is additively 
dissected into two quantities, D\ and D^,. Figure Q] illustrates a schematic picture of our 
framework for measuring the temporal turnover discussed above. 

The proposed turnover measure has the following properties: 


1. —oo < D < oo; 


2 . 

3. 

4. 


D = 0 <= D i = 0 and D -2 = 0 <;=> p(t) = p(u) and A (t) 
D\ < 0 (c/. Kullback and Leiblerl ( 1951 )): 

T>2 > 0 A (it) > A (t), T>2 < 0 A (it) < A (t). 


X (n); 


Note that Property 2 holds when the instantaneous change rate (Equation [T]) of the every 
species in the community is equal to zero, oti(t) = 0,2 = 1 , 2 ,..., s. 
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Figure 1: A schematic figure of the framework for measuring temporal turnover with four 
species (Spl, Sp2, Sp3 and Sp4) in a community at times t, t + 1,..., t + 5. For clarity, the 
x-axis represents observed-abundance, but note the fact that zero observed-abundance does 
not necessarily mean zero expected-abundance in the framework (for details see Section [3]). 
The proposed turnover measure D = D\ + D 2 (Equation [4|) quantifies the turnover of the 
species community due to two key ecological processes: change in community composition 
(D\) and change in total abundance (D 2 ). For times t to t + 3 — a new species Sp4 comes 
in at t + 1. both Sp3 and Sp4 drop out at t + 2 and only Sp3 comes back at t + 3 — the 
turnover is quantified by D\ as there are changes in community composition but not in total 
abundance (D 2 = 0). Between times t + 3 and t + 4, the situation is opposite: there are 
changes in total abundance but not in community composition. Hence, the turnover of this 
period is quantified by D 2 as D\ = 0. Both community composition and total abundance 
change between t + 4 and t + 5, the turnover of which is quantified by D\ and D 2 both. 
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The implications of Equation @ are worth consideration. First, it suggests that the 
temporal turnover is additively decomposed into two parts: the first term (D i) related to the 
amount of change in community composition, and the second term () being dependent only 
on the amount of change in community size sensu its abundance. This fact highlights two 
important aspects in evaluating the turnover of species community: 1) change in community 
composition and 2) change in total abundance. 

Second, it is interesting that the Kullback-Leibler divergence ( Kullback and Leibler, 19511 ) 
is subsequently derived from the definition (Equation [Sj) as D \. Kullback-Leibler divergence 
quantifies the difference between two relative abundance distributions of the community at 
different time occasions, p(t ) and p(u). Thus, the first term, D 1 , is interpreted as the part 
evaluating the compositional change of the community. 

The basis of our framework, equation ([2]) , highlights a key feature of temporal turnover 
the arrow of time. If A (t + h ) and A(f) are swapped, the instantaneous change rate cq(t) 
takes a different form. In other words, the interpretation of the temporal turnover measure 
is asymmetrical to time — a feature that contrasts with spatial turnover. Although the 
proposed turnover measure, D, could be applied to spatial contexts, by assigning t and u to 
different locations, it is important to note that its interpretation will differ from the temporal 
one. This is because of the asymmetrical structure inherent in temporal turnover whereas 
symmetricity is usually assumed when studying spatial turnover. 


3 Specifying the expected abundance A j(t) 

We introduced the idea of the expected abundance earlier, but have not specified it yet. 
Since we can only observe species abundances and never observe the expected abundance in 
ecological investigations, we need to determine the expected abundance from the observations. 
Here, we denote the observed abundance of the i-th species at time t as lower case, 

and treat it as a realisation of a random variable IVj(f), upper case, which follows the Poisson 
distribution with the time varying mean parameter Aj(f). We adopt the convention of using 
lower and upper cases for non random and random variables respectively. The expected 
abundance is then specified as the expectation (or the mean) of the random variable as 
A i(t) = E [jVj(t)]. A range of approaches are available for estimating the expected abundance, 
A i(t), from observations, {n,(t)}; the choice is largely dependent on the researcher and data 
characteristics, such as the number of observations. We give here a short summary of common 
approaches, widely used in population and community ecology. 

Many methods of estimating the expected abundance rely on the maximum likelihood 
principle. For instance, recalling the assumption we made, an approach widely used in the 
ecological literature is to use the observation itself as A i(t) = Here, the hat sign means 

its estimate. Using the observation as the estimate of the expected abundance can be regarded 
as the maximum likelihood estimate (MLE) when we have only one observation obtained for 
estimating the expected abundance at time t. Note that this may violate some properties of 
the turnover measure discussed in Section [2l when rii{t) = 0 is involved. Another choice is 
taking the average of the observations by assuming the expected abundance to be constant, 
A i(t) = A i, over the observation period, t = 1,2the MLE of which choice is then 
given as A, = J2t=i rii(t)/t n . 

An alternative approach is to model expected abundance, Aj(f), in relation to relevant 
environment factor(s). Although this approach requires more data than those outlined 
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above, a range of r egression type models are a vailable. Examples include generalised 
linear models (GLMs; McCullagh and Nelder . 19891 ) and generalised additive models (GAMs; 
Hastie and Tibshirani . 1990l ). that delineate the relationship b etween species abunda n ce an d 


environmental factors as well as other species’ abundances ( Kedem and Fokianos . 20021 ). 
This approach also entails specifying the link function and the distribution function that 
the abundance, Nj(t ), follows, wit h com m on ex amples including Poisson, Negative Binomial 
and others, see Zuur et al. ( 20091 ). Hilbe ( 2011 ) for details. This type of modelling approach 
offers more detailed insights in terms of interpreting community change as is discussed in the 
next section. We stress that the choice of estimation approach will depend on the extent 
to which information is available for study, and that m(t ) = 0 does therefore not necessary 
mean A*(t) = 0 unless the observation is used as the estimate. Many other approaches are, 
of course, applicable and not limited to those discussed above. For example, determining 


the tr a jectory of abundance, Nj(t), by (stochastic) differential equations, as did lives et al 
l200.il ) iMutshinda et all (hood l201ll ). is also a sensible approach. 


4 Identifying the influential drivers 


An advantage of modelling species expected abundance, A*(f), is that it becomes possible to 
identify drivers that influence the turnover measure, D. Equation Q reveals how this can 
be done. We illustrate this using a common model class, GLMs with the log link function, 
such as Poisson and Negative Binomial distributions, as an example. 

Suppose that the expected abundance of the i-th species, A*(f), is modelled by a GLM 
with the log link as 

m 

log(Aj(t)) = (5) 

3 =i 

where {xj(t)} are environmental variables, and {/%} are the parameters to be estimated. 
From Equations Q and Q, we have the turnover measure, D, described in an additive form 
as 

m s 

D(t : u) = ~ x j( t ))Pi( t )- ( 6 ) 

j =l *= l 

Since the contribution of the species and the environment factors are all additive, putting 

m 

di(t : u ) = ^2Pij(xj(u ) - Xj(t))pi(t ) and 

j=i 

S 

dj ( t:u ) = fa ( w ) ~ x 3 (*)» 

i =1 


the contribution ratio of the i-species and of the j-th environment variable to the turnover 
measure, D, is respectively defined as 


n(t : u) 


rj(t : u ) 


| di(t : u )| 
ZU\di(t : u)\' 
I dj(t : u)\ 

Z?= 


(7) 

( 8 ) 
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These quantities show what proportion each factor contributes to the absolute amount of 
the turnover. Whether this type of additive decomposition is available is largely dependent 
on how the species expected abundance, Aj(f), has been modelled. Although our example 
makes use of a GLM, as long as the right-hand-side of model ([5]) is additive in terms of the 
environment variables, as it is in GAMs, the additive decomposition (Equation[6|) is available, 
and consequently the contribution ratios, Equations 0 and (f8j). can be obtained. A key 
fact is that the turnover measure D can be additive in terms of species i and environment 
factor j. given an appropriate model form and a link function. A similar idea can be found 
in calculating the contribution of individu al species b ased on t he Bray-Curtis index (for 
example, simper of an R package ‘vegan’, Oksanen et ah . 20141 ) although, calculating the 
contribution of environment variables may not be straightforward since the Bray-Curtis index 
is, in general, non additive with respect to the environment variables. 


5 Simulation study 

To illustrate how our turnover measure works, and the kind of information that can be 
obtained, we perform a simulation study, and compare t he results with tho se prod uced using 
two p opular turnover metrics: Jaccard ( Jaccard . 1 190 ill and Bray-Curtis ( Bray and Curtisl. 
19571 ). We simulate an ecological community within whic h each species h old inter-species 
relationships keeping the zero-sum community condition ( HubbeM . 200 ll ). apart from 
imposed abrupt change to represent a disturbance as per Dornelas ( 2Q10h . 


an 


5.1 Simulated data 

In the simulation study, we assume an ecological community of 30 species (s = 30), within 
which species abundances at time t, N(t) = (W (t), N?,(t ), ..., N x (t)) T , are dependent on 
each other in a sense of zero-sum community ( Hubbelll . 20011 ). The abundances distribute as 
a multinomial distribution, IV(f) rsj Mn (n(t),p(t)), and the covariance between any pair of 
species, i and i !, is given as Cov [Ni(t), !Vj/(i)] = —n(t)pi(t)pi' (t) for i ^ i'. We generate series 
of abundance values (200 time steps for each of the 30 species), {N(t) = n(t) : 1 < t < 200}, 
for the whole community with an artificial change taking a place at t = 100 as 


Pr {N(t) = n(t)-,n,p) = 


Mn(n,p) t < 100, 
Mn(n',p') t > 100, 


where n and n' are the total abundance of the community and p and p' are the relative species 
abundances. In other words, the simulated community is a zero-sum community whose total 
size is respectively n and n' before and after the time point t = 100. The initial values of 
t he re lative species abundances, p , are obtained using the program provided by Dornelas! 
( 2Q10l ). We consider three types of changes in community size at time t: stable (n = n'), 
an increase (n < n') and a decrease (n > n'). For the compositional change, we consider 
two cases: no-alteration (p = p') and alteration (p ^ p') in the relative species abundances. 
The scenarios to be considered are specified by a combination of those states, but we have 
omitted the trivial case, stable and no-alteration (n = n' and p = p '), so that there are five 


scenarios to be examined. See Appendix C for the detailed explanation. 
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5.2 Calculating turnover metrics 

For each of our five simulated community time series, {rii(t) : 1 < t < 200 ,i = 
1,2,..., 30; n, n',p,p'} (n = 200, n' = 180, 200,230), the three turnover metrics: our measure 
D, Jaccard and Bray-Curtis indices, are calculated in two different ways. One is based on a 
GLM and the other is based on the simulated raw data, A i(t) = ruit). 

First, we fit a GLM with a Poisson distribution as 


log (Ai(i)) = On + 0 i2 I{t > 100 ), 


( 9 ) 


where On and 0*2 are the parameters to be estimated by maximising the likelihood. The 
fitted model (Equation 0 is exact re-parametrisation of the simulation data to capture the 
abrupt change introduced at time t = 100. Once the expected abundance is estimated, 
A i(t), the estimated relative abundance pi(t) can be calculated as pi(t) = Aj(f)/A(f), where 
A (t) = Ylt= iAj(i). The turnover metrics: our turnover measure, Jaccard and Bray-Curtis 
are then calculated. However, these indices are originally defined on observations, {rii(t)}, so 
we introduce a model-based version of Jaccard (J\) and Bray-Curtis ( BC \) indices for this 


study. See Appendix B 


Second, we also calculated the three turnover metrics based on the simulated data, in 
other words the observation rij(f) itself, A i(t) = rii(t). In doing so, we adopt the convention 
log(0/pi(u))0 = 0 for D\ of Equation (|3|) when the z-th species is newly observed in the 
community at time u, and also we omit the case where the i-th species becomes absent at 
time u, \og(pi(t)/0)pi(t) = oo, although this treatment may conceal the cases D\ = — oo, it 
is still useful for illustration purposes. Alternatively, an arbitrary very small value, e ~ 0, 
could be used for A i(t) when rii(t) = 0 to avoid numerical difficulties. 


5.3 Result 

For illustration purposes, here we show the result (Figure [2]) of the scenario, a decrease in 
the total abundance (n = 200 > n' = 180) and an alteration (p ^ p') in the relative species 
abundances. See Figures lAllIA4I in Appendix C for other scenarios. Since we are interested 
in documenting cumulative changes, all turnover metrics are calculated against the first point 
of the time series (i.e. the baseline). To reduce sensitivity to the exact nature of the baseline 
chosen, taking an average over some time points could also be a possible option. 

In this simulation, the Jaccard index of dissimilarity, calculated on the raw simulated 
data (1 — J) and on the model (1 — J\), reveals no turnover o f the com munit y str ucture. 
In contrast, the Bray-Curtis index, BC (original metric as in iBrav and Curtis! (119571 )) and 
Appendix B) does detect a shift in turnover. This contrast reflects the fact that 


as m 


BC X _ 

species richness of the simulated community has remained relatively constant throughout, 
the dominant species are observed constantly but the low abundance species are observed 
occasionally, whereas the abundances of the species within the community shifted at the 
perturbation point, t = 100. This analysis also highlights differences between the metrics 
based on the observations and those based on the model. The model-based ones tend to be 
robust against the variation due to observations, a feature that could be useful for detecting 
this kind of abrupt change. 

Our new turnover measure offers additional insights relative to the Jaccard and 
Bray-Curtis indices discussed above. D\ (the green line) quantifies the change in the 
community composition. With D\ the absolute divergence from zero indicates change has 
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Figure 2: Simulation results of the scenario with a decrease (n = 200 > n' = 180) and 
alteration (p ^ p'). All turnover measures are calculated against the first time point. The 
top row shows the relative abundances before (a) and after (b) the artificial change introduced 
at t = 100; (c) the simulated respective abundance series over 200 time steps; (d) the new 
turnover measure, D (red), D\ (green) representing the change in the community composition 
and D 2 (blue) representing the change in the community abundance. The solid lines are 
calculated based on the model (Equation [9j) , and the dashed lines are based on the simulated 
data; (e) the outcome of Jaccard index, 1 — J (dashed line, simulated data based) and 1 — J\ 
(solid line, model based); (f) the outcome of Bray-Curtis index, BC (dashed line, simulated 
data based) and BC\ (solid line, model based). 
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occurred in the relative abundance distribution of the community. For example, the addition 
of an invasive species to a community, leading to the local extinction of a native species 
would be detected by this metric. The change in community abundance is indicated by D 2 
(the blue line). When D 2 is positive the community abundance has increased, and when 
it is negative community abundance has decreased. For those scenarios involving decreased 
community abundance, n > n' (Figures HT1 and [7ul in Appendix C), D 2 is actually negative 
suggesting the community abundance decreased, whereas the Bray—Curtis index struggles to 
show such a difference, whether the community abundance is increased or decreased, because 
of it accounts for absolute differences. Once again the model-based measure of turnover, D 
and Z? 2 , provide a clear indication of the perturbation. A mass mortality event associated 
with a chemical spill, for example, would be detected by this metric. 


6 Application 

6.1 Data 


We further explore the performance of our new turnover measure using an exceptionally 
complete estuarine community time series. Sampling took place at Hinkley Point ‘B’ power 
station on the southern bank of the Bristol Channel in Somerset, UK (51°14 , 14.05 // N, 
3°8 / 49.71 // W). Monthly quantitative sampling of a fish community commenced in January 


2012) (Henderson and Holmes. 

1991, Henderson and Seabv, 

Henderson and Bird, 

2010. Henderson et all 201 ll). Ambient envir 


2005, Henders on. 2007, 


for details of the survey and its methodology. 


6.2 Model 


Following IShimadzu et al. ( 20131 ). we fit the same GAMs to the 45 c ore species that are 
consistently present in the assemblage ( Maeurran and Henderson . 2003h . The GAM for each 
7-th abundant species is 


log(Aj(t)) = P i0 + /ii(Year) + f i2 { Tide, height) + / i3 (Water. temp) + / i4 (Month), 

where the /5*o is a constant and f tJ (•) is a smoothing spline function. For the remaining 36 
species that occur infrequently we fit a constant model as the simplest model, 

l°g(A i(t)) = P i0 - 

Once the expected abundance is estimated, A i(t), the estimated relative abundance pi(t) can 
then be calculated as pi(t) = Aj(f)/A(f), where A(f) = )>^i=i A,(f). 


6.3 Result 

For illustration purposes, we present the turnover result for the month August over the period 
(Figured]). The top panel (Figure Eli) shows the extent to which the fish community changes 
each year relative to the first observation, August 1981. The turnover measure D (the red 
line) shows a reasonably stable fluctuation a little below the zero level over the period, and 
also suggests cyclical variation, since the late 1990s. The two components of the turnover 
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Figure 3: Illustrative results of the turnover analysis on the fish community at Hinkley 
Point in August over the period (1981 - 2012). Each year is compared against the first 
observation, August 1981. (a) The turnover measure ( D : the red line) and its components, 
the composition change {D\: the green line) and the community size change ( D 2 : the blue 
line); (b) The contribution ratios of each species in the community. Top five species with 
high contribution ratio are listed in the legend, and other colour/shading represent(s) the 
remaining 76 species; (c) the contribution ratios of each environmental factors: year (blue), 
tide height (green) and water temperature (yellow). 
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measure, one related to the change in community composition D\ (the green line) and the 
other related to the change in community abundance D 2 (the blue line), tell a more detailed 
story. The green line (D 1 ) shows the changes in the relative abundance distribution of the 
community, quantified as the departure from zero. The continuous decrease indicates that 
the species composition of the community has been gradually changing since 1981. On the 
other hand, the blue line ( D 2 ) resembles the cyclic variation of the red line ( D ) but with an 
increase over the period. This indicates that the community size has been getting slightly 
larger, relative to the 1981 level, since the late 1990s. 

In Figure [3}r, the contribution ratios of some species resemble the cyclical variation in D 
and D 2 . The middle panel illustrates each species’ contribution ratio defined by Equation ([7]) . 
Each coloured area represents the proportion of contribution to the change of the turnover 
measure D made by the species. A closer look reveals at least five species to be influential 
in the community: Trisopterus minutus, Pomatoschistus minutus, Merlangius merlangius, 
Solea solea and Sprattus sprattus (see the legend in Figure [3jo). The key role of these five 
species is reflected in their large contribution ratio (Equation [7]), inducing the departure 
of the turnover measure D from the zero level (Figure [3jr). The large area in Figure Or 
assigned to Trisopterus minutus means that it has been particularly influential over the 
period. Interestingly, Pomatoschistus minutus and Merlangius merlangius appear to be offset 
— as one increases in contribution ratio the other decreases, and their cyclical variations seem 
to correspond to the turnover measure D and D 2 in the late 1990s. Sprattus sprattus , a known 
dominant species in the community, has a largely constant contribution since the late 1980s. 

The bottom panel (Figure |3fc) shows the contribution ratio of the three environmental 
factors: year (blue), tide height (green) and water temperature (yellow). The influence of 
each factor varies through time. The peaks in the year effect match with those of D and D 1 , 
implying the fact that this community has been under the strong influence by some particular 
speices whose abundance fluctuates in a relatively regular cycle. 

7 Relationship with ct-diversity measures 

The framework that we have proposed is conceptually linked with a-biodiversity measures. 
Measuring biodiversity using Shannon’s entropy is equivalent to examining community 
change rate, but between two particular occasions: the reference time and the time when 
the community composition has become uniform. Consider a case where the community 
consisting of s species at the reference time t and at a particular time u, ( u > t ) when the 
relative abundance distribution of the community is homogeneous, p(u) = 1/s. The turnover 
measure D (Equation [2]) is then reduced to 

S 

D(t : u) = - log(s) - 5>(t) log (jpi{t)) + log 

i= 1 

If the community abundance stays the same, A(f) = A (u), the turnover measure is further 
simplified as 

S 

D(t : u) = Di(t : u) = - log(s) - ^Pi(f) log (jPi{t)). 

i= 1 

This is exactly Shannon’s entropy with a constant shift, — log(s). In other words, the two 
terms on the right-hand-side are, respectively, a species richness measure and a heterogeneity 
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measure. Inter estin gly, th is links back to the introduction of Shannon’s entropy for measuring 
biodiversity by Margate! (1957, 1958h . Combining the above reasoning with our framework 
discretising differential equations, it is clear that Margalef’s idea of using Shannon’s entropy 
bridges theoretical and empirical work within biodiversity research. Theoretical population 
ecology has modelling of the change rate of populations as its heart. As we have explained 
above, empirical studies that calculate Shanon’s entropy are doing something that is 
mathematically equivalent to investigating the change rate, but using field data. 

There have been many attempts to achieve generality in the heterogeneity component 
of biodiversity, although by doing so t he theoretical link with the community change rate 


as mentioned earlier is lost. S 3 (Il973h is one of the first researchers to notice that 


community heterogeneity can be measured by changing the emphasis on species domin ance 
in a c ommunity. In terms of the turnover measure, D, a generalisation equivalent to Hill 
ljl973h can further be made on the change in community comp osition, D\. Since Hill’s 


number is the natural log scale of Renyi’s entropy (IRenvil 


y compo 

IB, 


Renyi’s divergence that 


encompasses Kullback-Leibler divergence as a special case may be used as a measure of the 
change in community composition, D\. There hav e also been a va ri ety of divergence measure s 
examined in information t heory, for example see Kawada (1987), Read and Cressie ( 19881 ). 
Cichocki and Amari ( 201C)h . 


/3-diversity is often discussed in the context of par titioning of divers i ty with other types 


of diversity measures, such as 7 - and ct- diversities (IWhittakerl . Il960l . lLandel. Il996l . I. lost! . 


_ __ 1 i _ _ —i l _. ^- - •> - 

20071 . IChao et all l2012l ). In this context iMarcon et al.l ( 2012l ) have derived a /3-diversity 


measure, based on S hann o n’s e ntropy and its decomposition, using weighted Kullback-Leibler 
divergence. Reeve et al. 163) have also studied diversity partitioning in the spatial context, 
and have discussed a /3-diversity measure that considers species similarity by taking the 
exponent of Kullback-Leibler divergence. However, we stress that our turnover measure, D, 
has been derived directly from the community dynamics. Moreover the derivation of D is 
independent of biodiversity partitioning. 


8 Bray—Curtis index and its decomposition 


Temporal /3-diversity (temporal turnover) represents the change in biodiversity between two 
occasions. The Bray-Curtis index is an abundance-based metric of community dissimilarity, 
initially developed to measure spatial change but now increasingly used in the temporal 
context. Although it is, to be precise, defined using observed abundances, not on the 
expected abundance, there is a direct link between the turnover measure we propose and 
the Bray-Curtis index (see Appendix B for the original definition). In fact, the Bray-Curtis 


index, BC, can be rewritten as a function of our turnover measure, D, as 


BC = 


n=ii*i 

2 + D 


( 10 ) 


where \d t \ is the absolute value of di = \og(Xi(u)/Xi(t))pi(t), and recall that the turnover 
measure D is the sum of these di s over the species i ( i = 1,2, ...,s) which make up 
the community, so that D = Yli=idi (Equation [3]). Figure [4^ illustrates this functional 
relationship /Equation HOD . The vertical striped area shows where the value of Bray-Curtis 
index (BC) can lie, given the value of the turnover measure ( D ), as the numerator of Equation 
m varies depending on the extent to which the species in the community increase or 
decrease. The black solid line represents the lower bound of Bray-Curtis index, the case 
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Figure 4: (a) Functional relationship of the turnover measure, D, and Bray-Curtis index, 
BC. The vertical striped area, {(D, BC) : Y2t= 1 |dj|/(2+-D)}, is the values of the Bray-Curtis 
index can take, given the value of the turnover measure. The boundary (the black solid line) 
represents BC = \D\/{2 + D ) as |D| < Ylt= i |; (b) The relationship between the balanced 
variation in abundance, dsc-bah and the change in the species composition, D\\ (c) The 
relationship between the abundance gradient, dsc-gra, and the the change in the community 
abundance, Z? 2 - 


where the abundances of all species have increased. The lower bound clearly shows that 
the Bray-Curtis index is asymmetrical against the community abundance change. In other 
words, the Bray-Curtis index tends to amplify the effect of decreases in community abundance 
com pared t o its increase . 


Baselga (l2010l . 2012 , 20 bf) stud ied a decomposition of /3-biodiversity metrics, such as 
Jaccard, Sprensen ( Sorensen! 19481 ) and Bray-Curtis indices, and showed that they can 
be additively decomposed into two components. In particular, the Bray-Curtis index can 
be dissec ted into balanced variation in abundance and the abundance gradient components 
(Baselga, 2013h . We note here that our decomposition (Equation (H) has parallels with but is 
not equivalent to these components. In fact, the balanced variation in abundance, dsc-bal-, 
and the abundance gradient, dsc-gra , are functions of the turnover measure, D = D\ + D 2 ■ 
The balanced variation in abundance is described as 


dBC-bal 


'D+e~ D2 , (D < 0) 
-D~, (D> 0) 


( 11 ) 


where D + = Yh=i (M“) “ Aj(t)) + /A(t) and D = Y,t=i(^i(. u ) ~ A*(t))_/A(t). Here (•)+ 
and (•)_ respectively denotes its positive or negative part; note that D = D + + D~ but 
clearly D\ 7 ^ D2 7 ^ D + or D~ . This means that the balanced variation in abundance can 
be interpreted as the proportion of abundance increased or decreased. Figure 2)3 illustrates 
that there is little relationship between the balanced variation in abundance, dsc-bah and 
the change in the species composition, D\. 
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The abundance gradient is also described as 


dBC—gra 


\D\ 


2 + D 
D 

2 + D 


qe 


q, 


-d 2 


(D< 0) 
(D> 0) 


( 12 ) 


where q = min(Aj(u), Xi(t))/X(t). This is proportional to Bray-Curtis index itself, and 
is its lower bound. Figure illustrates the relationship between the abundance gradient, 
dsc-gra, and the ch ange in th e tot al abundance, D2. For the detailed derivations of Equations 
Baselga ( 20131 ) paper, 


(fill) and 


m 

drel 


from 


Leeendrel (120141 1 has recently 


see 
reviewed 


Appendix A 


approaches to partitioning /3-diversity, and 
identifi ed two major classes: the Baselga family tha t we h ave discus sed above, and the Podani 


family (jPodani and Schmeral . l201ll . ICarvalho et al 


13) 


(2013) for the discussion on partitioning 


2013l l: see also iLegendre and De Cace^l 
/3-diversity. Since the Podani family is 


mathematically related to the the Baselga one , it can a lso be described as a function of our 


turnover D. The R package ‘BAT’ (Cardoso et, all 1201 4l ) can be used to partition /3-diversity. 


9 Conclusion 


We have developed a framework for measuring the temporal turnover of species communities 
based on community dynamics resulting from processes, such as local immigration, extinction 
and pop ulation grow th. Our framework decomposes turnover into two classical ecological 


rules ( Browni . 1981): change in the community composition associated with allocation 
rules and change in community abundance due to capacity rules. The new temporal 
turnover measure combines these two aspects in an additive manner. The formal study 
of the framework highlights the links with other widely used methods of measuring a- and 
/3-diversity. 

The performance of the new turnover measure was examined by a simulation study and the 
analysis of an estuarine fish community. A comparison with other common indices, namely the 
Jaccard and Bray-Curtis indices, demonstrated that our measure offers new insights into the 
ecological basis of community change. In addition the simulation study has clearly shown that 
the approach we have described in this paper can provide useful information in investigating 
the community change based on observational data. Furthermore, we have demonstrated that 
our framework offers an important practical advantage in allowing researchers to identify the 
species and/or environmental factors that play important roles in community change through 
time. 

Taken together, these features of this new approach suggest that it will be a useful 
framework for quantifying and interpreting temporal turnover i n ecologica l communities at 
a time when the natural world is facing unprecedented threats ( Butchart et al. . 2010 ). 
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Appendix A The derivation of Equations (TTTT) and (TTZl) 


Putting A i(u) = Xij and A i(t) 


Xj.k, we have eqn 1 


3 in iBaseleal (j2013h 


as 


A = y^min(Aj(^), Xj(t)) 

i =1 

s s 

B = y A i(u) - min(Ai(u), A;(f)) = ^(Aj(u) - A i(t))+ 

2—1 2—1 

S S 

C = y A j(t) - min(Aj(u), A j(t)) = - y~](A j(u) - A$(£))_ 

2=1 2=1 

where (■)+ and (•)_ respectively denotes its positive or negative part. The Bray-Curtis index 

J2i=i I K{u) — Aj(t)| _ B + C 
A i{u) + A i(t) 2 A + B + C 
min (B,C) \B — C\ A 

A + min(B, C) 2A + B + C A + min(B, C) 

dBC-bal + dBC-gra- 

The balanced variation for each case is respectively given as 


is then described as 

dBC = 


Case I: min(S, C) = B (B < C <=> D < 0) 

B _ Ei=i(M“) - Ai(t))+ 


dBC-bal ~ 


A + B 


\{u) 


= Ei=l(Aj(lt) - Aj(t)) + A(f) n + -Da 
A (i) A (it) ' ’ 


Case II: min (B, C) = C (C < B D > 0) 


dBC-bal ~ 


C 


£-=i (Xi(ti)-Ai(t))_ 


= -£>- 


A + C A(i) 

The abundance gradient for each case is respectively given as 
Case I: min(S, C) = B (B < C D < 0) 

\B — C\ A _ | A (it) - A(t)| 22i =1 min(\(u),\(t)) 


dBC—qra — 


2A + B + CA + B 


A (it) + A (t) A(it) 

IA (it) — A(2)| A (t) y =1 mm(Xi(u), Xi(t)) 


= \D\e 


X(t) A(it) A(it) + A(f) 

-d 2 _S_ 


2 + D' 


Case II: min(S, C) = C (C < B D > 0) 

B-C A 


dBC—qra — 


A(it) - A (t) £- = 1 min(Aj(u), A;(t)) 


2A + B + CA + C 


= L> 


A(u) + A(t) 

q 


X (t) 


2 + D' 
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Appendix B Model-based Jaccard and Bray—Curtis indices 


For the simulation study, we calc u late two common indices: Jaccard ( Jaccardl . Il90lh 
and Bray-Curtis ferav and Curtisl. Il957h indices for comparison. As these indices are 
originally defined on observations, (nj(t)}, we introduce a model-based version of Jaccard 
and Bray-Curtis indices as follow. 


Jaccard index (J) 

j = ELi-fM*) > 0LML > 0) 
ELi + n i( u ) > °) 

Model-based Jaccard index ( J \) 


T, 8 i=iP(Nj(t) > 0 )P(Nj(u) > 0) = ELi d - exp(-A,(t))) (1 - exp(-A j(u))) 
ELi + N i(u) >0) s - ELi exp(-A i(t) - Xi(u)) 


Bray-Curtis index (BC) 

BC = ELi I ni(t)-m(uy 

n(t) + n(u ) 

where n(t) = ELi 

Model-based Bray Curtis index ( BC \) 


BC X = 


ELi IMt) -\j(u)\ 

A (t) + A('u) 


where A (t) = ELi E(0- 

Podani et al. ( 2013l l note a link between Jaccard and Ruzicka ( Ruzicka . 1958l l indices, 


taking m(t ) = 1 for species presence and rii(t) = 0 for its absence. The same link also exists 
between Sprensen and Bray-Curtis indices, but not between those, Jaccard and Bray-Curtis 
indices that we have used. This type of continuity is not crucial in our study since the 
relative ordering is, whether Jaccard or Sprensen (say L ) index is used, retained because of 
their monotonic relationship, L = 2 J/(J + 1). 
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Appendix C Simulation results 


We have studied the five simulation scenarios with the three states in changes of community 
size: stable (n = n'), an increase (n < n') and a decrease (n > n')\ and of community 
composition: no-alteration (p = p') and alteration (p ^ p'), when an artificial change 
taking a place at time t = 100. The scenarios are specified by a combination of those 
states, but we have omitted the trivial case, neutral and no-alteration (n = n' and p = p'), 
so that there are five scenarios to be examined. The scenario 3 is presented in the main text. 


Scenario 

Abundance change 

Composition change 

1 

n > n' 

p = p' 

— 

n = n' 

p = p' 

2 

n < n' 

p = p' 

3 

n > n' 

P + P' 

4 

n = n' 

P + P' 

5 

n < n' 

P + P' 


For example, the interpretation of Scenario 1 here is that the total community size 
(abundance) decreases n to n' by an artificial change taking a place at time t = 100, but the 
community composition, the relative species abundances p, stays as the same (p = p'). 
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Figure Al: Simulation results of the scenario with a decrease (n = 200 > n' = 180) and 
no-alteration (p = p'). All turnover measures are calculated against the first time point. 
The top row shows the relative abundances before (a) and after (b) the artificial change 
introduced at t = 100; (c) the simulated respective abundance series over 200 time steps; (d) 
the new turnover measure, D (red), D\ (green) representing the change in the community 
composition and D 2 (blue) representing the change in the community abundance. The solid 
lines are calculated based on the model (Equation [9|) , and the dashed lines are based on the 
simulated data; (e) the outcome of Jaccard index, 1 — J (dashed line, simulated data based) 
and 1 — J\ (solid line, model based); (f) the outcome of Bray-Curtis index, BC (dashed line, 
simulated data based) and BC\ (solid line, model based). 
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Figure A2: Simulation results of the scenario with an increase (n = 200 < n! = 230) and 
no-alteration (jp = p'). All turnover measures are calculated against the first time point. 
The top row shows the relative abundances before (a) and after (b) the artificial change 
introduced at t = 100; (c) the simulated respective abundance series over 200 time steps; (d) 
the new turnover measure, D (red), D\ (green) representing the change in the community 
composition and D 2 (blue) representing the change in the community abundance. The solid 
lines are calculated based on the model (Equation [9|) , and the dashed lines are based on the 
simulated data; (e) the outcome of Jaccard index, 1 — J (dashed line, simulated data based) 
and 1 — J\ (solid line, model based); (f) the outcome of Bray-Curtis index, BC (dashed line, 
simulated data based) and BC\ (solid line, model based). 
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Figure A3: Simulation results of the scenario with being stable (n = n' = 200) and alteration 
(p ^ p'). All turnover measures are calculated against the first time point. The top row shows 
the relative abundances before (a) and after (b) the artificial change introduced at t = 100; (c) 
the simulated respective abundance series over 200 time steps; (d) the new turnover measure, 
D (red), D\ (green) representing the change in the community composition and D 2 (blue) 
representing the change in the community abundance. The solid lines are calculated based 
on the model (Equation [9|) , and the dashed lines are based on the simulated data; (e) the 
outcome of Jaccard index, 1 — J (dashed line, simulated data based) and 1 — J\ (solid line, 
model based); (f) the outcome of Bray-Curtis index, BC (dashed line, simulated data based) 
and BC\ (solid line, model based). 
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Figure A4: Simulation results of the scenario with an increase (n = 200 < n! = 230) and 
alteration (p ^ p'). All turnover measures are calculated against the first time point. The top 
row shows the relative abundances before (a) and after (b) the artificial change introduced 
at t = 100; (c) the simulated respective abundance series over 200 time steps; (d) the new 
turnover measure, D (red), D\ (green) representing the change in the community composition 
and T >2 (blue) representing the change in the community abundance. The solid lines are 
calculated based on the model (Equation [9j) , and the dashed lines are based on the simulated 
data; (e) the outcome of Jaccard index, 1 — J (dashed line, simulated data based) and 1 — J\ 
(solid line, model based); (f) the outcome of Bray-Curtis index, BC (dashed line, simulated 
data based) and BC\ (solid line, model based). 
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Appendix D R code 


We provide R code to calculate our turnover measure and to produce the figures presented 
in the main text. 

DD 

Calculates the turnover measure, D , D\ and D 2 against any reference time point (the default 
is the first time point). 

Usage 

DD (x) 

DD(x, ref.t = 1, zero.rm = FALSE) 

Arguments 

x 

ref. t 

zero.rm 


Values 

A data frame that contains the turnover measure D, D\ and D 2 in each column. 

Example 

results <- list() 

for(l in l:n) results[[i]] <- glm(y xl + x2, family="poisson") 
x <- do.call(cbind, lapply(results, predict, type="response")) 

DD (x) 

Code 

DD <- function(x, ref.t=l, zero.rm=FALSE){ 
lmb <- apply(x, 1, sum) 

D2 <- log(lmb/lmb[ref.t]) 

x.p <- t(apply(x, 1, function(z)z/sum(z))) 

Pt <- x.p [ref.t,] 
if(zero.rm==FALSE){ 

D1 <- -t(apply(x.p, 1, function(z)ifelse(Pt==0, 0, log(Pt/z)))) y o *°/ 0 Pt 
}else-[ 

D1 <- -t(apply(x.p, 1, function(z)ifelse(Pt==0Iz==0, 0, log(Pt/z)))) %*% Pt 

} 

D <- D1 + D2 

data.frame(D, Dl, D2) 

> 


a data frame (t, i ) that contains the estimated expected abundances, 
A i(t). 

reference time point, t, the turnover measure is calculated against. The 
default value is ref . t=l, the first time point. 

a logical value indicating whether pi(u) = 0 should be stripped before the 
computation proceeds. The default value is zero. rm=FALSE. To avoid 
any NaN produced when observations are used as the estimated expected 
abundances, try zero. rm=TRUE. 
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cr.cal 


Calculates the contribution ratios of the species ri(t,u) or the contribution ratios of the 
environment factors rj(t,u ) from a list object, each branch of which must contain a GLM or 
a GAM object. 

Usage 

cr.cal(x) 

cr.cal(x, id=l, ref.t=l) 

Arguments 
x 
id 


ref. t 

Values 

A data frame containing the cumulative contribution ratios that can be passed to the 
function cr.plot to draw a contribution ratio diagram. 

Example 

x <- list() 

for(l in l:n) x[[i]] <- glm(y xl + x2, family="poisson") 
cr.cal(x) 


a list contains a GLM or a GAM output. 

a numerical value indicating whether the contribution ratios of the 
species, ri(t,u), (id=l) or the contribution rations of the environmental 
factors, rj(t,u), (id=2) to be calculated. The default value is id=l. 
reference time point, t, the turnover measure is calculated against. The 
default value is ref. t=l, the first time point. 


Code 


cr.cal <- function(x, id=l, ref.t=l){ 

lmb <- do.call(cbind, lapply(x, predict, type="response")) 
lmb.p <- t(apply(lmb, 1, function(z)z/sum(z))) 

lmb <- lapply(x, function(z)as.matrix(predict(z, type="terms"))) 
lmb <- lapplyClmb, function(z)t(t(z)-z [ref.t,])) 

lmb <- array(unlist(lmb), dim=c(nrow(lmb[[1]]), ncol(lmb[[1]]), length(lmb))) 


if(id==l){ 

absDi <- abs(apply(apply(lmb, c(l,3), sum), 1, lmb.p [ref.t,])) # sp 

}else{ 

absDi <- t (abs (apply (lmb, c(l,2), " 0 / 0 * 0 / 0 ", lmb .p [ref . t,]))) # env 

} 

ri <- apply(absDi, 2, function(z)z/sum(z)) 
ri.cum <- t(apply(ri, 2, cumsum)) 

data.frame(ri.cum) 
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cr.plot 

Plots a contribution ratio diagram (eg. Figures [3h> and[3]:-). 

Usage 

cr.plot(Y) 

cr.plot(Y, x=NULL, col.pal=NULL, k=5, ...) 

Arguments 
Y 

x 

col.pal 
k 


Example 

Y <- cr.cal(x) 
names(Y) <- col.names 
cr.plot(Y) 


a data frame (t,i) or (t,j) that contains the contribution ratios of each 
species, ry(t, u), or each environmental factor, rj(t,u). 
the coordinates of points in the plot. 

the colour palette to be used. The length needs to be the same as 
ncol(Y). The default choice is rainbow. 

the number of species or environment factors to be listed in the legend 

of the figure. The default value is k=5. 

further arguments passed to or from other methods. 


Code 


cr.plot <- function(Y, x=NULL, col.pal=NULL, k=5, . ..){ 
k <- min(ncol(Y), k) 

if(is.null(col.pal)){col.pal <- rainbow(ncol(Y))} 

Y <- data.frame(0, Y) 

if(is.null(x)) x <- c(l:nrow(Y)) 

matplot(x, Y, type="n", ylab="Contribution Ratio", ...) 
for(j in 2:ncol(Y)){ 

polygon(c(x, rev(x)), c(Y[,j-1], rev(Y[,j])), col=col.pal[j-1], border=NA, 

> 

lgd.list <- apply(apply(Y, 1, diff, na.rm=T), 1, sum, na.rm=T) 
legend("topright", legend=names(rev(sort(lgd.list)))[l:k], 
col=col.pal[order(lgd.list, decreasing=T)][1:k], 
pch=15, pt.cex=2, bg="white", ...) 
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